home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / int.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  2.5 KB  |  80 lines  |  [TEXT/????]
  1. (load "/scmutils/math/vector.scm")
  2.  
  3. (define (integrate f initial-states dt n integrator)
  4.   (let ((step-map (integrator f dt)))
  5.     (let loop ((n n) (states initial-states))
  6.       (if (= n 0)
  7.           states
  8.           (loop (-1+ n) (cons (step-map states) states))))))
  9.  
  10. (define (forward-euler f dt)
  11.   (lambda (states)
  12.     (let ((state (car states)))
  13.       (add-vectors state
  14.                    (scalar*vector dt
  15.                                   (f state))))))
  16.  
  17. (define (backward-euler f dt)
  18.   (lambda (states)
  19.     (let ((state (car states)))
  20.       (let ((predicted
  21.               (add-vectors state
  22.                            (scalar*vector dt
  23.                                           (f state)))))
  24.         (add-vectors state
  25.                      (scalar*vector dt
  26.                                     (f predicted)))))))
  27.  
  28. (define (trapezoid f dt)
  29.   (lambda (states)
  30.     (let* ((state (car states))
  31.            (d (f state))
  32.            (predicted
  33.              (add-vectors state
  34.                           (scalar*vector dt d))))
  35.       (add-vectors state
  36.                    (scalar*vector (/ dt 2)
  37.                                   (add-vectors (f predicted) d))))))
  38.  
  39. (define (milne f dt)
  40.   (let ((predictor (trapezoid f dt)))
  41.     (lambda (states)
  42.       (let* ((sn (car states))
  43.              (sn-1 (cadr states))
  44.              (predicted (predictor states)))
  45.         (add-vectors sn-1
  46.                      (scalar*vector (/ dt 3)
  47.                                     (add-vectors (f predicted)
  48.                                  (scalar*vector 4 (f sn))
  49.                          (f sn-1))))))))
  50.  
  51. (define (circle state)
  52.   (vector (- (vector-ref state 1))
  53.           (vector-ref state 0)))
  54.  
  55. (define fe
  56.   (reverse (integrate circle (list (vector 1 0)) .1 200 forward-euler)))
  57.  
  58. (define be
  59.   (reverse (integrate circle (list (vector 1 0)) .1 200 backward-euler)))
  60.  
  61. (define tr
  62.   (reverse (integrate circle (list (vector 1 0)) .1 200 trapezoid)))
  63.  
  64. (define mr
  65.   (reverse (integrate circle 
  66.                 (integrate circle (list (vector 1 0)) .1 1 trapezoid)
  67.                       .1
  68.                       199
  69.                       milne)))
  70.  
  71. (define (show)
  72.   (text-graphics 9)
  73.   (draw-graph fe 'per-point connect-the-dots 'screen-box (left-half-box))
  74.   (draw-graph be 'per-point connect-the-dots 'screen-box (right-half-box))
  75.   (clear-box (left-half-box))
  76.   (draw-graph tr 'per-point connect-the-dots 'screen-box (left-half-box))
  77.   (clear-box (right-half-box))
  78.   (draw-graph mr 'per-point connect-the-dots 'screen-box (right-half-box))
  79.   (restore-video-mode))
  80.